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Abstract 

We link the spontaneous breakdown of chiral symmetry in Euclidean QCD to the col- 
lision of spectral shock waves in the vicinity of zero eigenvalue of Dirac operator. The 
mechanism, originating from complex Burger's-like equation for viscid, pressureless, 
one-dimensional flow of eigenvalues, is similar to recently observed weak-strong cou- 
pling phase transition in large N c Yang-Mills theory. The spectral viscosity is propor- 
tional to the inverse of the size of the random matrix that replaces the Dirac operator 
in the universal (ergodic) regime. We obtain the exact scaling function and critical 
exponents of the chiral phase transition for the averaged characteristic polynomial for 
N c > 3 QCD. We reinterpret our results in terms of known properties of chiral random 
matrix models and lattice data. 



1. Introduction 

It was a recently emphasized that the Burgers equation could be used to understand 
universal features of the weak to strong-coupling transition in two-dimensional Yang- 
Mills theory with a large number of colors N c [ 1, Si. This transition, first studied 
by Durhuus and Olesen 0], can be pictured, in the language of the Burgers equation, 
as resulting from the collision of two spectral shock waves at the closure of the gap. 
The emergence of the Burgers equation in such kind of problems, and related questions 
in random matrix theory, seems generic, and can be simply understood by exploiting 
Dyson's original idea of matrix random walks. Then, after a proper rescaling of time 
that separates the fast motion of the eigenvalues caused by their mutual repulsion, from 
the diffusion that results from the random walks of the matrix elements, one can easily 
show, by using standard tools of statistical mechanics, that the average resolvent obeys 
indeed a Burgers equation (or its simple generalizations) in the limit of large size of 
the matrices. Here we shall obtain the Burgers equation from an exact equation, akin 
to a diffusion equation, satisfied by the average characteristic polynomial. The average 
resolvent and average characteristic polynomial are simply related in the large N limit. 
However, the equation for the characteristic polynomial is exact for any matrix size, 
while no such equation seems to exist for the average resolvent. 
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Returning to two-dimensional Yang-Mills theory, we note that the average charac- 
teristic determinant for a Wilson loop of a given area was indeed shown to obey an exact 
diffusion equation for any number of colors N c (which indeed yields a Burgers equa- 
tion in the infinite number of colors limit). For finite number of colors, the flow of the 
eigenvalues is viscid, with a negative spectral viscosity v s — *w. It is the negative sign 
of the viscosity, or of the diffusion constant, that causes the rapid spectral oscillations 
at the closure of the gap. This picture was confirmed by extensive numerical studies of 
full Yang-Mills theory in three dimensions and recently, in four dimensions yLl£|]. 

In this letter, we demonstrate that a similar mechanism, governed by similar viscous 
spectral Burgers-like equation, is responsible for the universal spectral oscillations of 
the spectrum of the Dirac operator in QCD that accompany the spontaneous breakdown 
of the chiral symmetry BaLZI]. In the case of the Durhuus-Olesen transition, where the 
role of the time is played by the area of the Wilson loop, an explicit construction of 
a random matrix model was provided by Janik and Wieczorek [8i] , matrices attached 
to large loops being formed by multiplying random matrices representing small loops. 
Here, as already mentioned, following Dyson |9|], we add a fictitious time (somewhat 
analogous to Schwinger's proper time) in order to describe the diffusion of the random 
matrices in 4 dimensional Euclidean space. 



2. Random matrix theory and QCD 

A cornerstone for the microscopic understanding of the spontaneous breakdown of 
chiral symmetry is the Banks-Casher II 1 Oil relation 



l<M>|sZ = »4r, (i) 

where the quark condensate < qq > is an order parameter for chiral symmetry, p(0) 
is the averaged (over the gauge field configurations) level density of the Euclidean 
Dirac operator near the vanishing eigenvalue, and V4 = L 4 is the Euclidean volume. 
This relation shows that chiral symmetry breaking requires a strong accumulation of 
eigenvalues near zero, i.e. a level spacing A ~ l/L 4 , much larger than the level spacing 
A ~ l/L of a free system [11]. This accumulation of eigenvalues leads to universal 
properties, that are well captured by random matrix theory: For eigenvalues smaller 
than a characteristic energy scale, referred to as the Thouless scale £77,, the fluctuations 
of eigenvalues are described by chiral random matrix models respecting the global 
symmetries of the Dirac Hamiltonian. In QCD, the condition E®, D /A = F^L 2 » 1, 
where F„ is the pion decay constant, determines the regime of applicability of random 
matrix theory ifTHflll . 

In Euclidean QCD, all four Dirac matrices can be chosen to be anti-hermitian, 
hence the spectrum of the massless Dirac operator D = iy^id^ - igA^) is purely imagi- 
nary. The partition function, for a fixed topological sector, reads 

Z? c » = lY]m}Y](A> + m 2 f )\ (2) 

W I ' y 
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where the averaging is done with respect to gluonic configurations of a given topologi- 
cal charge v, ±iA^ are the eigenvalues of D, and nif is the mass of a quark with flavor /. 
Due to the chiral symmetry non-zero eigenvalues of D come in pairs, and the number 
of fermionic zero modes is related to the topological charge. In the chiral Gaussian 
random matrix model (hereafter ^GUE), corresponding to QCD with N c > 3, the role 
of the massless Dirac operator is played by a random matrix W(- —iD) of the following 
form: 

Here K is a rectangular M x N (M > N) matrix with complex entries, K/j = x,j + /y, ; , 
where Xij and yy are drawn from a Gaussian distribution. Note that W is hermitian, so 
that its eigenvalues k,-'s are real. The block-diagonal structure of d3) reflects the chiral 
symmetry of the Dirac operator: W anticommutes with the analogue of the Dirac matrix 
75, defined here as 75 = diag(ljy, -1m)- This implies in particular that the eigenvalues 
come in pairs of opposite values, (k, —k). By construction, W has in addition v = M—N 
zero eigenvalues. These mimic the zero modes of quarks propagating in gauge fields 
of non trivial topology. 

General spectral properties of the random matrices can be obtained form correlation 
functions containing both products and ratios of the characteristic polynomial [ 14] 

/ nf =1 z(w ; ) \ 

C(w 1 ,....wf,w 1 ,....w b ) = l (4) 
\Uj=iZ(Wj)l 

where < ... > denote averaging with respect to the ^-GUE measure, and the charac- 
teristic polynomial is Z(yv) = det(w - W). Directly related to (|4]i is the resolvent, 
R(z) = d w C(z; h0| w=z , whose imaginary part yields the spectral density. In the vicinity 
of zero, the microscopic (unfolded) resolvent predicts for QCD 115I1 

^ = x(h(x)K h (x) + I M K b ^) (5) 

where b — Nf + |v|, with Nf the number of quark flavors, and x = zV^L. The appearance 
of Bessel functions If,, Kt, in the correlation functions, such as in Eq. Q, is generic. 
They encode the universal behaviors that show up already in the simplest objects, the 
average characteristic polynomial and/or the average of the inverse of the characteristic 
polynomial. As the first result of the paper, we shall obtain an exact differential equa- 
tion for the averaged characteristic polynomial. This equation is akin to a diffusion 
equation, from which a Burgers equation can be derived from a simple transformation. 



3. Burgers equation and the characteristic polynomial for Dirac operator 

We assume now that the entries of the matrix K follow independent random walks. 
Let us denote by P(X, Y, t) the joint probability that the entries of K take the values 
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X = {xij} and Y = {yij} at time t. The random walks of the matrix elements translates 
into the following diffusion equation for P(X, Y, t): 

j M N 

d,P (X, Y, t) = - ^ + d 2 Y )P (X, Y t) . (6) 

(=1 j=l 

We shall be interested in this paper with the time evolution of the averaged charac- 
teristic polynomial 



Q v n {w, t) = <det[w - W[) = j dXdY P(X, Y, t) det[w - W(X, Y)], (7) 

where n — N + M, w is an arbitrary complex number, and dX = Ylij dx !; , and similarly 
for dF. In order to get the equation obeyed by Q„{w, t), we consider first the equation 
satisfied by the average characteristic polynomial of the associated N x N Wishart 
matrix K^K, M V N = (detjz - K*k\). This, we write as an integral over Grassmann 
variables r\iJj{. 



M v N (z,t) = J di/d^e' 2 '* J dXdYP(X,Y,t) exp^-ljK^Kr]), 



(8) 



where TjK T Kri = 2 ;t=1 2 y=1 ^(^y _ iyij)( x jk + iyjk)Vk- To derive the equation for M, 
we take a time derivative of the expression above. This acts on P{X, Y, f), which, using 
Eq. ©, we transform into derivatives with respect to x, ; and y,j. Then, we integrate 
by parts, and use standard Grassmann calculus to obtain (after a somewhat tedious but 
straightforward calculation) the following differential equation 

d t M v N (z, t) = -zd 2 z M v N {z, t) - (y + 1) d z M v N {z, t). (9) 

This equation is valid for any and M, and arbitrary initial conditions. Note that for 
the trivial initial condition Kfj(t = 0) = 0, its solution is given by time dependent 
associated Laguerre polynomial iflfill . 

From the equation for M v N (z, t), Eq. (0 above, one easily obtains the equation for 
Q v N (w, f) = w v M^iz = w 2 , t). This reads 

d,Q v „(w, t) = -i#G>, f) - ^d w Q v n (w, t) + ^e;( W , f). (10) 

It will be also useful to consider the equation for the Cole-Hopf transform of Q v n {w, t), 
j% = ifl w ln(Q^(w, t)). This object identifies with the average resolvent in the large n 
limit. After a rescaling of the time, r = Mt, one gets from Eq. ( TTOb 



■^d T f: + r n d w f: = ~ 

n In 



\n) w J 



(11) 



where we have separated on the left hand side the terms that survive the large n limit, 
and on the right hand side the terms that are explicitly suppressed by powers of 1/n. 
Note the crucial role played by the rescaling of time in arriving at this equation. The 
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motivation behind this rescaling is that the diffusion associated with the random walks 
is taking place over a time scale that is larger, typically by a factor n, than the time 
scale corresponding to the local rearrangements of the eigenvalues due to their mutual 
repulsion 09j,ll6|]. After rescaling, the diffusion terms are dwarfed by a factor 1 /n, and 
the large n dynamics is dominated by repulsion. The last term, of order 1 /n 2 finds it 
origin in the kinematical zero modes present when v + 0. 



4. Large n limit 

We consider now the limit n — > oo, with v constant. We set g(w , t) = lim^oo /„ v (w, t). 
Eq. dTTb reduces then to the inviscid Burgers equation Jl, 17, 18], independent of v: 



d T g(w, r) + g(w, T)d w g(w, t) = 0. (12) 

It can be solved using complex characteristics. We choose the system to be initially in 
a chiral symmetric state, with eigenvalues localized at +a.The characteristic lines are 
given by w = £, + Tgo(^), with 

g0 ( w ) = g ( w , T = 0) = U— + -L-) = — " (13) 
2 \ w — a w + aj w z — a L 

The solution g(w, t) is constant along the characteristics, meaning g(w, t) = go [£(w, t)] . 
By eliminating one obtains an implicit equation for g: 

tV -2-rwg 2 + (w 2 -a 2 + r)g -w = 0. (14) 

This equation can be solved by elementary means, and well known results recovered. 
In fact, the change of variables a> = w/ a/t, d = aj s/t, and G s ,(cj) = -v/t^(w\/t, t) 
transforms this equation into G 3 st - 2a>G 2 t + (to 2 - d 2 + l)G sr = cj, an equation for a 
time-independent resolvent G sr that has been obtained in this context using different 
techniques 11191 ud, 12111 . In previous studies, the parameter d was introduced as the 
"deterministic" part of the chiral matrix (with K in Eq. (O replaced by K + d, d being 
fixed and K random), in order to control the approach to the chiral transition. The 
aforementioned change of variables renders transparent the dynamics captured by the 
Burgers equation: For small time, i.e. t < a 2 , the spectral density remains localized 
in humps centered around the values ±a. As time reaches the critical value t c = a 2 
(corresponding to d — 1 in the static approach), the two domains of the spectrum 
merge at the origin, which we picture as the collision of two spectral shock waves. A 
finite "condensate" then develops, and chiral symmetry is spontaneously broken. 

We now recover these features by studying the singularities of the characteristics, 
and the behavior of the solution in the vicinity of these singularities. This brief dis- 
cussion will pave the way for the scaling analysis to be performed in the next section. 
Singularities appear when characteristics start to cross (appearance of a spectral shock 
wave). This occurs for values of £ that obey the equation, 

^ =0=1+t^(&), (15) 
04r f=£ 
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that is 

r(a 2 + ^ 2 )^(f c -a 2 f. (16) 

The values of £ c correspond to the edges of the spectrum when r < t c = a 2 . At the 
critical time t c = a 2 , where the two humps of the spectrum start to merge, the equation 
for f c admits a double solution at £. = 0, which splits into two purely imaginary and 
opposite solutions when t > a 2 . 

In order to study the behavior of g at the edge of the spectrum, we expand go in the 
vicinity of a singular point 

go® = go®) + (? - + \(£ - tcfg'o®) + \<£ - + • ■ • (17) 

When t < a 2 , g' ®) = -I It and, using go = (w - £)/t, we get 

w-w c = ^-&) 2 ^'(&) + ... (18) 

One can then easily invert the relation between w and and get (for the rightmost 
edge) 



^•^te^' <19) 



so that 



g(w, t) = go© - go(ft) + (? - &)g (&) = go©) + - J^tt V^ 3 ^, (20) 

t y^gd®) 

which exhibits the familiar square root behavior of the spectrum near its (right) edge. 

For t - a 2 ,a similar analysis taking into account that go(^c) = = g'^®), so that 
the cubic term must be kept, yields (for w c = 0) 

*<w.t) = • (2D 

For time r > a 2 , we have g(w = 0, r) = go(£c(w = 0, r)) = - ^~ T , which is 
imaginary and hence directly proportional to the spectral density p(0). 

The behavior of g(w) at the edge of the spectrum determines the average eigenvalue 
spacing in the limit of large matrices. A singularity ~ \w - w c \" yields a level spacing 
~ rT s with 6 = 1/(1 + a). We have therefore 6 = 2/3 for r < a 2 , 5 = 3/4 for r = a 2 
and 6 = 1 for t > a 2 . We shall exploit these properties in the next section. 



5. Critical properties of the characteristic polynomial 

In this section we carry out a scaling analysis of the average characteristic polyno- 
mial Q(w, t), or its Cole-Hopf transform f(w, t), in the vicinity of the singular points. 
To that aim, we set 

w = w c + n s s, f N -» go®) + n^xis, r), (22) 

with y = 1 -5, and s and^(i, r) remain finite as n — > oo. 
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5.0.1. Airy edge 

Let us focus first on the left edge of the positive part of the spectrum, for r < a 2 , 
where we expect the solution to be of the form (see Eq. (l20l >) 



/",>,T)*Sofe)+Sofe). 



2(w- w c ) 
rg'^c) ■ 



(23) 



This exhibits a square root singularity (5 - 2/3), and in line with the discussion above 
we set w — w c + rC^ s and / w v — > go(£c) + n~~x( s > T )- (Notice that go(fr) is a function 
of t only, go(^ c ) - 5 T w c .) Substituting this ansatz into Eq. ( TTTI i. and keeping only the 
dominant terms as n — > oo, we get (with go(£c) = 3 T go(fc)) 



go 



(&) + \d] X = = «9 P (2/tf (£) +* 2 + a*) ■ 



(24) 



This equation is easily integrated, x 2 + d s x + 2sgo(£ c ) + = 0, with m(t) an arbitrary 
function of r. Defining x( s > T ) = 3, In </>(i, t) and shifting s = S - u(t) / (2go(^ c )) (with 
0(s, t) = ^(S, r)) one transforms the equation above into the equation for the Airy 



function, cni/' + 2go(£ e )/?i/r = 0. The looked for solution is therefore 



0(s, t) = Ai 



-(2g )Hs + —) 
2go 



(25) 



The arbitrary function u(t) can be determined by a careful matching of the asymptotic 
form as i — > oo of the solution to its large N limit. One then finds that u actually 
vanishes 

5.0.2. Bessel universality 

For r > a 2 , there is no singular behavior (5 = 1). So, we set w = n~ x s and j% — x- 
In the large n limit (at v constant) we obtain, following the same manipulations as 
above, the following partial differential equation: 



= d s [x 2 + d sX +~-^ 



(26) 



which integrates to^ 2 + d s x + f 
obtain: 



+ m(t) = 0. Then, setting xis, t) = d s \n(p(s, t) we 
s 2 d](f> + sd s cf> + cf> [s 2 u(t) - v 2 ] = 0, (27) 



whose solution is 



(s, t) = /„ [s V M ( T )] • 



(28) 



The determination of the arbitrary function m(t) proceeds as in the previous case, by 
matching the asymptotic x( s > t) ~ —i V«(t) with the large AT solution. This gives 
V«(t) = Vr - a 2 /T. We recover the scaling of the ratio of spectral densities discussed 
for instance in 11221 l23l 12411 . 
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5.0.3. Bessoid (axially symmetric Pearcey) universality 

Finally, we move to the case of r — t c — a 2 , which will lead to the second new 
result of this paper. As we have shown above, at r = a 2 , the two pre-shocks collide. 
Before the collision, these pre-shocks are accompanied by oscillations of the Airy type. 
Our purpose now is to describe the modification of the pattern of oscillations for t 
close to t c . To perform this analysis, it is most convenient to start from the diffusion 
equation obeyed by the characteristic polynomial, i.e., Eq. (TTOb which, after changing t 
into t = Mt, we rewrite as 

^ G >' T) = -W dlQV " M ~ 4^' G « (W ' T) + 4^ e> ' T) ' (29) 

This equation is to be solved with the initial condition Q v n (w,r — 0) = w v (w 2 - a 2 ) 
(where the function of w, (w 2 - a 2 ) N is defined with a cut between -a and a). It can 
be verified by a direct calculation (that does not require the explicit calculation of the 
integral below), that 



e>,r) = Cr 



'f 

Jo 



y exp M — 



/v (^) (/ _ fl Ydy, (30) 



where C = (-1) V+I ^ 2 2M, is a solution with the proper initial condition. (For v = 0, 
this solution agrees with a known solution B25I1 .) The y-integral runs over a half-line 
that starts at the origin and goes to infinity, making a constant angle </>,. = arg(y) with 
the real axis, with — n < arg(y) < n. For the integral to be convergent, we require | < 
\<f>y\ < 2ze. The modified Bessel function I v (x) had the following asymptotic expansion 



I v (x) a —p=e x , valid for I arg(x)l < ? (see 12611 ; here x - 2Myu ). This is useful 

\2ttx 1 T 

in particular to verify the initial behavior. Indeed, as r — > 0, one may estimate the 
integral using the saddle point method. The saddle point equation yields y + w — 0, 
which fixes in particular arg(y) = arg(-w). This new condition for <p y , together with the 
convergence condition noted above, are easily seen to be compatible with the condition 
of validity of the asymptotic expansion of the Bessel function. In turn, these conditions 
limit the allowed arguments of w to | < | arg(w)[ < ^f. 

The integral representation (|30t of the characteristic polynomial allows us to study 
the vicinity of the critical point. We note that the saddle point equation reads (in the 
large n limit) 



2y 2w 

T T 



y 



J 7 

y- — a z 



= 0. 



(31) 



Identifying y = — £ we recognize the equation for the characteristic lines. This indi- 
cates how the large n dynamics is coded in this integral. We shall focus more specif- 



ically at the critical point, w 
ln(a 2 - y 2 ) * ln(a 2 ) - : -" 



Q,y = 0,t = a . In this regime, we may expand 



^4, and obtain 



a 2N C t- y exp 



exp(/0 v )co 



y v+I exp 



Ny 4 



Ny z My 2 



dy, (32) 
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with C = C(-l) N . To capture the critical behavior also as a function of time, we set 
r = a 2 + 8. Then Eq. (f32t becomes 



a 2(N-l) n , 



(Mw 2 \ r a m>> 00 . [ Ny 4 f My 2 e] (2Myw\ , , oos 

This expression suggests the following change of variables that will ensure a smooth 
large n limit: 

_ 1 _3 _ 1 2 

y = n *au, w — n *raq, t — n ^ra 8. (34) 

We then define 

B v {q,8) = lim C'~ l a- n n^(-\y +l Q v n \tf*raq, a 2 + -«"ira 2 6»|, (35) 

N,M— >°° \ 1 J 

where the limit is taken with v constant. Finally 

M v+1 expl--M 4 - -u 2 e\l v (qu)du. (36) 

This exact scaling function for the characteristic polynomial is the second important 
result of this paper. It has a form very similar to the Pearcey function, which gives the 
asymptotic behavior of the characteristic polynomial in the case when gap closes for 
GUE ll27tl or in the case of unitary diffusion on the circle |Ht] 



P(q,i)= [ dyexp(-y 4 -ty 2 +qy). (37) 

Here q is the rescaled angle representing the position of the eigenvalue on the unitary 
circle and t parameterizes the fluctuations around the critical area. The critical indices 
that determine the scaling with n in Eq. (l34l are identical to those in the Pearcey inte- 
gral, but the form of the integral is different. The reason is the chiral symmetry, which 
imposes an additional polar symmetry of the spectrum, trading the exponential function 
of q in the Pearcey integral (l37t for a Bessel function in Eq. 



6. Conclusions 

In this letter we have obtained an exact differential equation for the average char- 
acteristic polynomial of a chiral random matrix, and its Cole-Hopf transform. For the 
latter, the equation takes the form of a generalized viscid Burgers equation, where the 
viscosity is proportional to the inverse of the size of the matrix, but with a negative 
sign. This allowed us to provide a complete description of the full critical behavior 
of averaged characteristic polynomials in chiral QCD, based on a single equation. In 
particular, we considered the case of chiral Gaussian Unitary Ensembles and we have 
identified the exact universal scaling function (Bessoid B v (q, 8)) in the vicinity of the 
chiral critical point for the average characteristic polynomial. We did not analyze in 
this letter the properties of the average of the inverse characteristic polynomial, but we 
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have checked that if fulfills similar equations, albeit with initial conditions singular at 
w = 0, alike in the cases of unitary and GUE diffusions. 

Since F„ scales like V^o more and more eigenvalues of the Dirac operator fall 
into the universal window when the number of colors tends to infinity, the volume of 
the lattice being kept finite. This suggests, that the present study is relevant also for 
analyzing the spontaneous breakdown of chiral symmetry at finite volume and large N c 
QCD, which was observed, and explained by Neuberger and Narayanan 112811 . For small 
lattice sizes, chiral symmetry is unbroken, while at some critical scale L c a condensate 
is formed. The same authors [29] have also observed, that the Af c -dependence of the 
level spacing closest to zero goes from 1 /N c in the broken chiral symmetry phase to 
1 /N^ 3 in the symmetric (gapped) phase. At L — L c , the critical scaling changes to 
1 /N 3 ^ 4 behavior and the condensate vanishes at the critical size L c as y/L - L c . These 
results are in agreement with our analysis. 

The critical universal scaling function (Bessoid) for large N c Dirac operator resem- 
bles closely the critical universal scaling function (Pearcey's cuspoid) for the weak to 
strong coupling transition in Yang-Mills theory at large N c . It would be interesting to 
study numerically both transitions simultaneously (at least in some simple model like 



12811) to see the interplay between the cuspoid and the Bessoid, or, in other words, the 
relation between the critical size L c for chiral symmetry breakdown and the critical 
area ~ L? c for the weak to strong coupling transition in Yang-Mills theory. Of particu- 
lar interest would be to measure on the lattice the microscopic spectral density exactly 
at the point of the transition. As far as we know, the microscopic spectral density at 
criticality was constructed explicitly only for b — Nf + \v\ = case 13011 , and was never 
checked by the lattice simulation. Taking into account the still ongoing discussion on 
the nature of chiral phase transition, its relation to confinement and Anderson local- 
ization lf3lll32Tl . lattice verification of analytic predictions for microscopic densities at 
the critical point may be a powerful tool to shed more light on this aspect of strong 
interactions. 
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